Excess mortality: combining three countries

Data

expected_deaths_monthly <- read_rds("data/expected_deaths_monthly.Rds") %>% 
  mutate(Pandemic = case_when(
    (Year >= 1889 & Year <= 1892) ~ 1890,
    (Year >= 1917 & Year <= 1920) ~ 1918,
    (Year >= 1956 & Year <= 1959) ~ 1957,
    (Year >= 2019 & Year <= 2020) ~ 2020
  ))

Figure 1 - diffs

Monthly

Yearly

Figure 2 - four pandemics

Absolute

focus_years <- c(1890, 1918, 1957, 2020)

# plot_years <- c(seq(1890 - 1, 1890 + 2),
#                 seq(1918 - 1, 1918 + 2),
#                 seq(1957 - 1, 1957 + 2),
#                 seq(2020 - 1, 2020 + 2))

for(country in unique(expected_deaths_monthly$Country)) { 
  
  print(paste0("Plotting ", country))
  
  max <- expected_deaths_monthly %>% 
    filter(Country == country) %>% 
    filter(!is.na(Pandemic)) %>% 
    summarise(deaths = max(deaths),
              lpi_flu_hc = max(lpi_flu_hc),
              upi_flu_hc = max(upi_flu_hc)) %>% 
    rowwise() %>% 
    mutate(max = max(c(deaths, lpi_flu_hc, upi_flu_hc)))
  
  print(paste0("Max for x is ", max$max))
  
  for(year in focus_years) {
    
    print(paste0("    Year ", year))
    
    plot_data <- expected_deaths_monthly %>% 
      filter(Country == country) %>% 
      filter(Pandemic == year)
    
    name <- paste(country, year, sep = "_")
    
    if (year == 2020) {
      my_title = paste0(country, ": ", 
                        as.character(year - 1), "-", as.character(year))
    } else {
      my_title = paste0(country, ": ", 
                        as.character(year - 1), "-", as.character(year + 2))
    }
    
    plot <- plot_data %>% 
      ggplot(aes(x = Date)) +
      geom_rect(xmin = as.Date(paste0(year, "-01-01")),
                xmax = as.Date(paste0(year, "-12-31")),
                ymin = -Inf, ymax = Inf, 
                fill = "#ffffbf", alpha = 0.01, inherit.aes = FALSE) +
      geom_ribbon(aes(ymin = lpi_flu_hc, ymax = upi_flu_hc), alpha = 0.1) + 
      geom_line(aes(y = fit_flu_hc), colour = "grey50") +
      # geom_line(aes(y = median_flu_year), colour = "#2c7bb6") + 
      geom_line(aes(y = deaths), color = "#d7191c", size = 0.5) + 
      scale_x_date(limits = c(min(plot_data$Date), max(plot_data$Date)),
                   date_breaks = "6 months", date_labels = "%b") +
      xlab("") + 
      # ylab("Observed, 10y median & predicted deaths") +
      ylab("Observed & predicted deaths") +
      # labs(caption = "Excludes extreme years") +
      ggtitle(my_title) +
      scale_y_continuous(limits = c(0, max$max),
                         labels = label_number(accuracy = 0.1, suffix = "K", 
                                               scale = 1e-4)) +
      theme_bw() +
      theme(plot.title = element_text(size = 10),
            axis.text.x = element_text(size = 8),
            axis.title.y = element_text(size = 9),
            axis.text.y = element_text(size = 8))
    
    assign(name, plot)
  }
}
[1] "Plotting Switzerland"
[1] "Max for x is 10766"
[1] "    Year 1890"
[1] "    Year 1918"
[1] "    Year 1957"
[1] "    Year 2020"
[1] "Plotting Sweden"
[1] "Max for x is 17278"
[1] "    Year 1890"
[1] "    Year 1918"
[1] "    Year 1957"
[1] "    Year 2020"
[1] "Plotting Spain"
[1] "Max for x is 163422"
[1] "    Year 1890"
[1] "    Year 1918"
[1] "    Year 1957"
[1] "    Year 2020"
plot_spacer() + Spain_1918 + Spain_1957 + Spain_2020 +  
  Sweden_1890 + Sweden_1918 + Sweden_1957 + Sweden_2020 +
  Switzerland_1890 + Switzerland_1918 + Switzerland_1957 + Switzerland_2020 +
  plot_layout(widths = rep(c(2, 2, 2, 1), 3)) +
  plot_layout(ncol = 4) 

ggsave("paper/four_pndemics.png", dpi = 300, width = 297, height = 210, units = "mm")

Percent

plot_years <- c(seq(1890 - 1, 1890 + 2),
                seq(1918 - 1, 1918 + 2),
                seq(1957 - 1, 1957 + 2),
                seq(2020 - 1, 2020 + 2))

plot_data <- expected_deaths_monthly %>% 
  filter(Year %in% plot_years) %>% 
  mutate(Percent = ((deaths - fit_flu_hc) / fit_flu_hc),
         Difference = ifelse(deaths > fit_flu_hc, "More", "Less")) %>% 
  group_by(Country, Pandemic) %>% 
  mutate(time = row_number()) %>% 
  ungroup()

plot_data %>% 
  ggplot(aes(x = time)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_vline(xintercept = 13, colour = "grey80") +
  geom_vline(xintercept = 24, colour = "grey80") +
  geom_rect(xmin = as.Date(paste0(year, "-01-01")),
            xmax = as.Date(paste0(year, "-12-31")),
            ymin = -Inf, ymax = Inf,
            fill = "#ffffbf", alpha = 0.01, inherit.aes = FALSE) +
  geom_col(aes(y = Percent, fill = Difference), size = 0.5) +
  facet_grid(vars(Country), vars(Pandemic), scales = "free_x") + 
  scale_x_continuous(limits = c(min(plot_data$time)-1, max(plot_data$time)+1),
                     breaks = c(1, 13, 24, 36, 48), 
                     labels = c("-12m", "Start", "End", "+12m", "+24m")) +
  xlab("") + 
  ylab("% above predicted deaths") +
  ggtitle("Excess deaths as % of observed") +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(values = c("#67A9CF", "#EF8A62")) +  
  guides(fill = guide_legend(reverse = TRUE)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(), 
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(), 
        panel.grid.minor.y = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "grey40")) 

ggsave("paper/four_pndemics_perc.png", dpi = 300, width = 297, height = 210, units = "mm")

Table 1

pandemic_year <- c("1890","1918","1957","2020")

table_1 <- expected_deaths_monthly %>%
  select(Country,Year, deaths, pop, fit_flu_hc) %>%
  filter(Year %in% pandemic_year)%>%
  group_by(Year, Country) %>%
  summarise(Expected_Mortality = round(sum(fit_flu_hc), 0),
            Deaths_num = sum(deaths),
            Cum_excess_death = round(sum(deaths)-sum(fit_flu_hc), 0),
            Percent_excess_death = round((Cum_excess_death/Expected_Mortality)*100, 1))%>%
  ungroup()

openxlsx::write.xlsx(table_1,"paper/table1.xlsx", row.names = FALSE)
Year Country Expected_Mortality Deaths_num Cum_excess_death Percent_excess_death
1890 Sweden 73077 81824 8747 12.0
1890 Switzerland 58397 61805 3408 5.8
1918 Spain 451779 695758 243979 54.0
1918 Sweden 78642 104591 25949 33.0
1918 Switzerland 50473 75034 24561 48.7
1957 Spain 282888 293502 10614 3.8
1957 Sweden 70064 73132 3068 4.4
1957 Switzerland 52280 51066 -1214 -2.3
2020 Spain 425484 494604 69120 16.2
2020 Sweden 88975 95411 6436 7.2
2020 Switzerland 67893 75971 8078 11.9

Appendix

## Crude
expected_deaths_monthly %>% 
  ggplot(aes(x = Date, y = deaths)) +
  geom_col() +
  geom_rect(xmin = as.Date(paste0(1890, "-01-01")),
            xmax = as.Date(paste0(1890, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  geom_rect(xmin = as.Date(paste0(1918, "-01-01")),
            xmax = as.Date(paste0(1918, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  geom_rect(xmin = as.Date(paste0(1957, "-01-01")),
            xmax = as.Date(paste0(1957, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  geom_rect(xmin = as.Date(paste0(2020, "-01-01")),
            xmax = as.Date(paste0(2020, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  facet_grid(vars(Country), scales = "free_y") + 
  scale_y_continuous(labels = label_number(accuracy = 1L, suffix = "K", scale = 1e-4)) +
  theme_minimal() +
  xlab("") + ylab("Deaths") +
  ggtitle("Monthly number of deaths")

Per population

expected_deaths_monthly %>% 
  ggplot(aes(x = Date, y = (deaths/pop)*10000, colour = Country)) +
  geom_line() +
  theme_minimal() +
  xlab("") + ylab("Deaths / 10,000 population") +
  ggtitle("Monthly number of deaths")

Country panel

Switzerland

Alt Viz

Stacked

Diffs 2